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Recently, we have shown how a colored-noise Langevin equation can be used in the context of 
molecular dynamics as a tool to obtain dynamical trajectories whose properties are tailored to display 
desired sampling features. In the present paper, after having reviewed some analytical results for 
the stochastic differential equations forming the basis of our approach, we describe in detail the 
implementation of the generalized Langevin equation thermostat and the fitting procedure used 
to obtain optimal parameters. We discuss in detail the simulation of nuclear quantum effects, and 
demonstrate that, by carefully choosing parameters, one can successfully model strongly anharmonic 
solids such as neon. For the reader's convenience, a library of thermostat parameters and some 
demonstrative code can be downloaded from an on-line repository. 
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Stochastic differential equations (SDE) have been used 
to model the time evolution of processes characterized by 
random behavior, in fields as diverse as physics and eco- 
nomics. In particular, the Langevin equation (LE) has 
been regularly applied in the study of Brownian motion, 
and has been used extensively in molecular dynamics 
(MD) computer simulations as a convenient and efficient 
tool to obtain trajectories which sample the constant- 
temperature, canonical ensembleii^. 

In its original form, the Langevin equation is based 
on the assumption of instantaneous system-bath interac- 
tions, which correspond to the values of the random force 
being uncorrelated at different times. A non-Markovian, 
generalized version of the LE arises in the context of 
Mori-Zwanzig theorj«2il. In this theory, if one consid- 
ers a harmonic system coupled with a harmonic bath, 
it is possible to integrate out the degrees of freedom of 
the bath. This leaves one with a linear stochastic equa- 
tion where both the friction and the noise have a finite 
memory. The conventional LE is recovered in the limit of 
a clear separation between the characteristic time-scale 
of the system's dynamics and that of the system-bath 
interaction. 

This class of non-Markovian SDEs has been extensively 
used to model the dynamics of open systems interacting 
with a physically- relevant bath (see e.g. Refsi^"— ). In- 
stead, our recent works^iS have used colored(correlated)- 
noise SDEs as a device to sample efficiently statistical 
distributions in molecular-dynamics (MD) simulations. 
These works aimed to show how a stochastic thermo- 
stat suitable for Car-Parrinello-like dynamics^ could be 
constructed, and to include nuclear quantum effects in a 
large class of problems at a fraction of the cost of path- 
integrals calculations^. In these applications the real dy- 



namics is lost, and one focuses only on the efficient cal- 
culation of static ensemble averages. 

In this Paper we discuss the practical implementation 
of the generalized Langevin equation (GLE) thermostat 
that we used in the two cases mentioned above. We 
also provide the reader with the analytical and numeri- 
cal tools needed to construct SDEs tailored to their own 
sampling needs. Throughout we take advantage of the 
dimensional reduction scheme, which allows to exploit 
the equivalence between a non-Markovian dynamics and 
a Markovian dynamics in higher dimensionality. In doing 
this, we supplement the physical coordinates with addi- 
tional degrees of freedom^, whose equations of motion 
are taken as linear, so as to simplify the formalism and 
analytical derivations. 

In the Appendices we recall some of the properties of 
multidimensional stochastic processe o^'^"" — , which are 
useful to our discussion, and present a short comparison 
of the GLE thermostat and the widely used massive Nose- 
Hoover chainai^"— . A simple FORTRAN90 code imple- 
menting our method to the dynamics of an harmonic os- 
cillator and a library of optimized thermostat parameters 
can be downloaded from an on-line repositoryi^. 



I. GENERALIZED LANGEVIN THERMOSTAT 

A. Markovian and non-Markovian formulations 

The Langevin equation for a particle with position q 
and momentum p, subject to a potential V{q)^ can be 
written as 



q =p 

P = -V'{q) 



QjppP 



bppi{t). 



(1) 



'Reprinted with permission from J. Chcm. Theory Comput., 2010, 
6 (4), pp 11701180. Copyright 2010 American Chemical Society. 



where S_{t) represent an uncorrelated, Gaussian- 
distributed random force with unitary variance and zero 
mean =^ 0, (?(t)^(0)) = 5{t)]. Here and what fol- 
lows we use mass-scaled coordinates. Furthermore, for 
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consistency, the friction coefficient (usually denoted by 
7) is here given the symbol a^p, while bpp is the intensity 
of the random force. In this notation, the fluctuation- 
dissipation theorem (FDT) reads b^^ = 2appkBT. If this 
relation holds, the dynamics generated by Eq. ([T]) will 
sample the canonical ensemble at temperature T^il^. 

As explained in the Introduction, in order to bypass the 
complexity of dealing with a non-Markovian formulation 
directly, we supplement the system with n additional de- 
grees of freedom s = {s^} which are linearly coupled to 
the physical momentum and between themselves. The 
resulting SDE can be cast into the compact form 



q =p 



-V'{q) 




Clrin \ ( P 

A 



'pp 



hp B ' I ^ 



(2) 



Here, ^ is a vector of n + l uncorrelated Gaussian random 
numbers, with (6(00(0)) = StjS{t). Clearly, Eq. ^ 
is recovered when n = 0. For an harmonic potential 
^{q) = ^^^0^7 Eqs. ([2]) are linear, and an Ornstcin- 
Uhlenbeck process is recovered whose time propagation 
can be evaluated analytically. In the non-linear case one 
can use the Trotter-decomposition to split the dynamics 
into a linear part, which evolves the {p, s) momenta, and a 
non-linear part, which evolves the Hamilton cquationsi^. 
This is facilitated by the fact that the dynamics of (p, s) 
alone is linear, and its exact finite-time propagator can 
be analytically evaluated (see Subsection II E[) . 

Here and in the rest of the paper, we adopt the same 
notation introduced in Rcfj^ to distinguish between ma- 
trices acting on the full state vector x = {q,p, s)"^ or on 
parts of it, as illustrated below: 
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The Markovian dynamical equations ^ are equivalent 
to a non-Markovian process for the physical variables 
only. This is best seen by first considering the evolu- 
tion of the {p, s) variables in the free-particle analogue 
of Eqs. The additional degrees of freedom s can be 
integrated away, and one is left with (see Refi^ and Ap- 
pendix \^ 



{p, s) {p, s) ) , the drift matrix Ap and the diffusion ma- 



trix Bp is given by^ 



ApCp 



^p-^p 



BpB 



(6) 



In Appendix 1X1 we show that setting Cp = ksT is suf- 
ficient to satisfy the FDT. In this case, Eq. ([5]) fixes Bp 
once Ap is given. FDT also implies that the colored- 
noise autocorrelation function H{t) = (C(t)C(O)) is equal 
to kBTK{t), whereas the more complex relation between 
K{t) and H{t), valid in the general case, is reported in 
Eq. (1X51) . 

Since there is no explicit coupling between the position 
q and the additional momenta s, one can check that ex- 
actly the same dimensional reduction can be performed 
in the case of an arbitrary potential coupling p and q, and 
that Eqs. ^ correspond to the non-Markovian process 



dV /■* 

p=-—-j K{t~s)p{s)ds + C{t). 



(7) 



In the memory kernel ([5]), A can be chosen to be a general 
real matrix, and can have complex eigenvalues, provided 
they have a positive real part. This results in a K{t) that 
is a linear combination of exponentially damped oscilla- 
tions. Therefore, a vast class of non-Markovian dynamics 
can be represented by Markovian equations such as ([2]) . 



B. Exact solution in the harmonic limit 

The thermostats typically used in MD simulations have 
a few parameters, that are chosen by trial and error. A 
thermostat based on Eqs. ([2]) depends on a much larger 
number of parameters, and hence the fitting procedure is 
more complex. It is therefore important to find ways to 
compute a priori analytical estimates so as to guide the 
tuning of the thermostat. 

To this end, we examine the harmonic oscillator, which 
is commonly used to model physical and chemical sys- 
tems. By choosing V{q) = ■^ij-i^q^ the force term in ([2]) 
becomes linear, and the dynamics of x = (q,p, s)'^ is the 



OU process x = — AqpX-|-Bgp^. In Eqs. ([2]) the s degrees 
of freedom are coupled to the momentum only. There- 
fore, most of the additional entries in A^p and Bgp are 
zero, and the equations for x read 



P 



K{t-s)p{s)ds + C{t) 



(4) 



where the memory kernel K(t) is related to the elements 
of Ap by 



K{t) ^ 2appS{t) 



T - 

apG 



\t\A- 



(5) 



Based on the fact that the the free-particle dynamics 
of (p, s) is an OU process, one also finds than the re- 
lationship between the static covariance matrix Cp = 







-1 



cj^ app Sip 



a„ 




(8) 



The exact finite-time propagator for Eqs. ([8]) can be 
computed, and so it is possible to obtain any ensem- 
ble average or time-correlation function analytically. Of 
course, one is most interested in the expectation values 
of the physical variables q and p. In particular, one can 
obtain the fluctuations ^g^) and (p^) and correlation 
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functions of the form ^(7^(t)g^(0)), which can be used 
to measure the coupling between the thermostat and the 
system. The resulting expressions are simple to evaluate 
but lengthy, and we refer the reader to Appendix [B] for 
their explicit form. 

One can envisage, using the estimates computed for an 
oscillator of frequency a;, to predict and hence optimize 
the response of a normal mode of a similar frequency in 
the system being studied. Furthermore, thanks to the 
properties of Eq. ([5|), one docs not need to perform a 
normal-modes analysis to turn this idea into a practi- 
cal method. Consider indeed a perfect harmonic crystal, 
and apply an independent instance of the GLE thermo- 
stat to the three cartesian coordinates of each atom. It 
is easy to see that, since Eq. is linear, and contains 
Gaussian noise, the thermostatted equations of motion 
are invariant under any orthogonal transformation of the 
coordinates. Therefore, the resulting dynamics can be 
described on the basis of the normal modes just as in or- 
dinary Hamiltonian lattice dynamics. As a consequence, 
each phonon will respond independently as a 1-D oscilla- 
tor with its own characteristic frequency. Thus, to tune 
the GLE thermostat, one only needs the analytical results 
in the onc-dimcnsional case, evaluated as a function of uj. 
The parameters can then be optimized for a number of 
different purposes, based solely on minimal information 
on the vibrational spectrum of the system under inves- 
tigation, without any knowledge of the phonons eigen- 
modes. 

The invariance properties of the GLE thermostat lead 
to additional advantages. For instance, we can contrast 
its behavior with that of Nose- Hoover (NH) chains, based 
on equations which are quadratic in p (see Appendix [C|) . 
As a consequence of the nonlinearity, the efficiency of an 
NH chains thermostat for a multidimensional oscillator 
depends on the orientation of the eigenmodes relative to 
the cartesian axes, an artefact which is absent in our case. 

Having set the background, we now turn to the de- 
scription of the various applications of Eqs. 



tively) : 



C. Efficient canonical sampling 



We first discuss the design of a GLE which can op- 
timally sample phase space. In this case, the target 
stationary distribution is the canonical ensemble, so the 
equations of motion need to satisfy the detailed-balance 
condition. Still, there is a great deal of freedom available 
in the choice of the autocorrelation kernel or, equiva- 
lently, in the choice of Ap and Bp matrices. These free 
parameters can be used to optimize the sampling effi- 
ciency. To this end, we must first define an appropriate 
merit function. Standard choices are the autocorrelation 
times of the potential and total energy {V and H respec- 



Tv 



TH 



{{V{t)-{V)){V{Q)-{V)))dt 

(9) 

{{H{t)-{H)){H{Q)-{H)))dt. 



In the harmonic case, these can be readily computed in 
terms of correlation times of and (see Appendix IB|) . 
and will depend on Ap and the oscillator's frequency 
Lo. For example, one easily finds that in the white-noise 
limit, with no additional degrees of freedom as in Eq. p]). 



TH (w) 



4^' 



TV (w) 



2a 



VP 



Ujpp 

2^' 



(10) 



Both response times are constant in the high-frequency 
limit, and increase quadratically in the low-frequency ex- 
treme of the spectrum. For a given frequency one can 
choose Opp so as to minimize the correlation time - thus 



enhancing sampling. It should be noted that Eqs. ([TO]) 
contain a "trivial" dependence on w, as one expects that 
sampling a normal mode would require at least a time 
of the order of its vibrational period. One can thus de- 
fine a renormalized k{uj) = [r(ti;)w] ^ as a measure of 
the efficiency of the coupling. In the white-noise case, 
K = 1 for the optimally-coupled frequency {ujh ~ app/2 
and iJv = o,pp: respectively), and decreases linearly for 
lower and higher values of w. 

While this result in itself provides a guide to choose 
a good value of the friction coefficient in conventional 
(white-noise) Langevin dynamics, we can enhance the 
value of k(w) over a broader frequency range, by using a 
colored-noise SDE. If we want to obtain canonical sam- 
pling, the FDT has to hold, so that Cp = ksT. We 
therefore consider the entries of Ap as the only indepen- 
dent parameters, since Bp is then determined by Eq. ([6|). 

In practice, we set up a fitting procedure, in which 
wc choose a set of frequencies Wi, distributed over a 
broad range {uJmin, '-^max) ■ For an initial guess for the 
thermostat matrix Ap we compute k(w) for each of 
these frequencies. We then vary Ap, so as to optimize 
mini K^uji), and aim at a sampling efficiency on the range 
(wmm, Wmoa;) which is as high and frequency- independent 
as possible. We will discuss this fitting procedure in more 
detail in Section |TT1 

In Figure [T] we compare the optimized k(ll)) for differ- 
ent frequency ranges and number of additional degrees of 
freedom. We find empirically that k(w) = 1 is the best 
result which can be attained, and that nearly-optimal 
efficiency can be reached over a very broad range of fre- 
quencies. This constant efficiency decreases slightly as 
the fitted range is extended, regardless of the number 
n of Si employed. For a given frequency range, however, 
increasing n has the effect of making the response flatter. 

Clearly this scheme will work optimally in harmonic 
or quasi-harmonic systems, and anharmonicity will in- 
troduce deviations from the predicted behavior. In the 
extreme case of diffusive systems such as liquids, one has 
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FIG. 1: Sampling efficiency as estimated from Eq. ([9} for an 
harmonic oscillator, plotted as a function of the frequency u;. 
The k{lo) curve for a white-noise Langevin thermostat opti- 
mized for ij = 1 [black, dotted lines, Ea. (|10|) ] is contrasted 
with those for a set of optimized GLE thermostats. The pan- 
els, from bottom to top, contain the results fitted respectively 
over a frequency range spanning two, four and six orders of 
magnitudes around uj = 1. Dark, continuous lines correspond 
to matrices with n — 4, and dashed, lighter lines to n = 2. 
The GLE curves correspond to the sets of parameters |kv_4- 2 1 
kv_2-2, kv_4-4, kv_2-2, kv_4-6, kv_2-6| which can be down- 
loaded from an on-line repositorjsi^ 

to ask the question of how much diffusion wiW be affected 
by the thermostat, especially since in an overdamped 
LE equation the diffusive modes are considerably slowed 
down (see e.g. Refi^). To estimate the impact of the 
thermostat on the diffusion, we define the free-particle 
diffusion coefficient D* a.s that calculated switching off 
the physical forces. Its value when a GLE thermostat is 
used is 

mD* 1 f°° , , , , 

ksT (p^) Jo (11) 

where we have assumed the FDT to hold. In practical 
cases, if an estimate of the unthermostated (intrinsic) 
diffusion coefficient D is available, one should choose the 
matrix Ap in such a way that D* ^ D, so that the ther- 
mostat will not behave as an additional bottleneck for 
diffusion. Equation has the interesting consequence 
that D* can be enhanced either by reducing the over- 
all strength of the noise, as in white-noise LE, but also 
by carefully balancing the terms in the denominator of 
Eq. (HI]). 

We have found empirically that for an Ap matrix 
fitted to harmonic modes over the frequency range 
i'-^mim ^max) , thc diffusion Coefficient computed by (|lip 




T2 < 



FIG. 2: A cartoon representing a two-thermostat setup, which 
we take as the simplest example of a stochastic process violat- 
ing the fluctuation-dissipation theorem. If the relaxation time 
versus frequency curves for the two thermostats are different, 
a steady-state will be reached in which normal modes corre- 
sponding to different frequencies will equilibrate at different 
effective temperatures. 



is D* Ki ksT /iOmin- Thls latter expression gives a useful 
recipe for choosing the minimal frequency to be consid- 
ered when fitting a GLE thermostat for a system whose 
diffusion coefficient can be roughly estimated. 



D. Frequency-dependent thermostatting 

Thc ability to control thc strength of the thermostat- 
system coupling as a function of the frequency - demon- 
strated above - points quite naturally at more sophis- 
ticated applications. For instance, one can apply two 
thermostats with distinct target temperatures and dif- 
ferent efficiencies k{ijj) (see Figure [5]). Obviously, such 
a simulation is not an equilibrium one, since energy is 
systematically injected in some modes and removed from 
others, but leads to a steady state that has useful prop- 
erties. Indeed, the normal modes will couple differently 
to the two thermostat, so that the effective temperature 
of each mode can be controlled as a function of oj. This 
two-thermostats example is just an instance of a broader 
class of stochastic processes, for whom the FDT is vi- 
olated. In general, we can relax the assumption that 
Cp — ksT, and for a given drift matrix wc can choose a 
Bp which is suitable to our purpose. 

Returning to the harmonic oscillator case, one can 
solve exactly the dynamics for a given choice of Ap, Bp 
and frequency uj. The resulting dynamics is performed 
in the n + 2-dimensional space defined by the variables 
{q,p,s), according to Eq. For a compact notation, 

we used thc full matrices Aqp and Bgp. The full Cqp(a;), 
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which defines the stationary distribution in the steady 
state, can be computed solving an equation analogous 
to ©: 

AqpCqp + CqpAqp = B^pB^p (12) 

One can tune the free parameters (Ap and Bp) so as to 
make the Cqqitj) and Cpp(a;) elements of the extended co- 
variance matrix as close as possible to the desired target 
functions (g^) {uj) and (p^) (w). 

In a previous paper— we applied this method to ob- 
tain (g^) (w) and (p^) (w) in agreement with the values 
appropriate for a quantum harmonic oscillator, and ob- 
tained a good approximation to the quantum-corrected 
structural properties in quasi-harmonic systems. Many 
other applications can be envisaged, which take advan- 
tage of frequency-dependent thermostatting. For in- 
stance, one could use this technique in accelerated sam- 
pling methods^ir— , which work by artificially heating the 
low-frequency modes, whilst keeping the other modes at 
the correct temperature. 

E. Implementation 

The implementation of a GLE thermostat in 
molecular-dynamics simulations is straightforward. 
Here, we consider the case of a velocity- Verlet integra- 
tor, which updates positions and momenta by a time 
step At, according to the scheme: 

p^p + V'{q)At/2 

q^q+pAt (13) 
p^p + V'{q)At/2. 

Eqs. (jl3p can be obtained using Trotter splitting in a 
Liouville operator formalism^l. In the same spirit one 
can introduce our GLE thermostat by performing two 
free-particle steps by At/2 on the (p, s) variablesi^: 

{p,s)^P[{p,s),At/2] 

p^p + V'{q)At/2 

q^q + pAt (14) 

p^p + V'{q)At/2. 
{p,s)^V[ip,s),At/2] 

At variance with thermostats based on second-order 
equations of motion such as Nose-Hoover, where a mul- 
tiple time-step approach is required to obtain accurate 
trajectorie a^^i^^ , this free-particle step can be performed 
without introducing additional sampling errors. The ex- 
act finite-time propagator for (p, s) reads: 

V [(p, s) , Atf = T{At) (p, s)^ + S{At)e (15) 

where ^ is a vector of n + 1 uncorrelated Gaussian num- 
bers, and the matrices T and S can be computed once, 
at the beginning of the simulation and for all degrees of 



freedoniiSiSi. The relations between T, S, Ap, Cp and 

At read 

T = e-'^*^^ SS^ = Cp - e-^*^-Cpe-^*<. 

It is worth pointing out that when FDT holds, the 
canonical distribution is invariant under the action 
of ([T51) . whatever the size of the time-step. A useful con- 
sequence of this property is that, in the rare cases where 
applying (jlSp introduces a significant overhead over the 
force calculation, the thermostat can be applied every 
m steps of dynamics, using a stride of m At. This will 
change the trajectory, but does not affect the accuracy 
of sampling. 

The velocity- Verlet algorithm (|13p introduces finite-At 
errors, whose effect needs to be monitored. In micro- 
canonical simulations, this is routinely done by checking 
conservation of the total energy H. Following the work 
of Bussi et al^ we introduce a conserved quantity H, 
which can be used to the same purpose: 

H = H AK, (16) 

i 

where AKi is the change in kinetic energy due to the 
action of the thermostat at the i-ih time-step, and the 
sum is extended over the past trajectory. In cases where 
the FDT holds, such as that described in Section II CI 
the drift of the effective energy quantitatively measures 
the violation of detailed balance induced by the velocity- 
Verlet step, similarly to RefsJ^i^. In the cases where 
the FDT does not hold, such as the frequency dependent 
thermostating described in Section [I Dl the conservation 
of this quantity just measures the accuracy of the inte- 
gration, similarly to Refs i^^'^° . 

II. FITTING OF COLORED-NOISE 
PARAMETERS 

A key feature of our approach resides in our ability 
to optimize the performance of the thermostat based on 
analytical estimates, making the method effectively pa- 
rameterless. Such optimization, however, is not trivial, 
even if computationally inexpensive. The relationship 
between Ap, Bp and the correlation properties of the re- 
sulting trajectory is highly nonlinear. Furthermore, we 
found empirically that many local minima exist which 
greatly hinder the optimization process. With these dif- 
ficulties in mind, we provide a downloadable library of 
fitted parametersii which can be adapted to most of 
the foreseeble applications, according to the prescriptions 
given in Section HID I Details about the fitting procedure 
are given in the following three subsections. 

A. Parameterization of GLE matrices 

A number of constraints must be enforced on the drift 
and diffusion matrices in order to guarantee that the re- 
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suiting SDE is well-behaved. It is therefore important 
to find a representation of the matrices such that during 
fitting these conditions arc automatically enforced, and 
that the parameters space is efficiently explored. A first 
condition, required to yield a memory kernel with expo- 
nential decay, is that all the eigenvalues of Ap must have 
positive real part. A second requirement is that the ker- 
nel K{uj) is positive for all real uj. This ensures that the 
stochastic process will be consistent with the second law 
of thermodynamics^. 

Finding the general conditions for Ap to satisfy this 
second constraint is not simple. However, we can state 
that a sufficient condition for K{uj) > is that Ap -I- A J 
is positive definite. For simplicity we shall assume such 
a positivity condition to hold, since we found empirically 
that this modest loss of generality docs not significantly 
affect the accuracy or the flexibility of the fit. Moreover, 
in the case of canonical sampling, Ap + Aj > is also 
required in order to obtain a real diffusion matrix, since 
BpBj ksT (Ap + Aj) according to Eq. ©. 

One would like to find a convenient parameterization, 
which enforces automatically these constraints. This is 



best done by writing Ap = A- 



(S) 



the sum of a 



symmetric and antisymmetric part. Since any orthogonal 
transform of the s degrees of freedom would not change 
the dynamics (see Appendix E^. one can assume without 

loss of generality that the A'^'^^ block in Ap'^' is diagonal 
(see Eq. ([3]) for the naming convention). Since in the 

(A) 

general case the antisymmetric Ap does not commute 

(S) (S) 

with Ap , we will assume it to be full, while Ap can 
be written in the form 



a 


ai 


a2 ■ 


■ an \ 


ai 


ai 
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a2 
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an 





• 


• an J 



(17) 



In order to enforce the positive-definiteness, one use an 



analytical Cholesky decomposition Ap 



(S) 



QpQp , with 
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• 


• dn J 



(18) 



and Ui = df, ai = diqi, and a = + J^ilf- Such a 
parameterization guarantees that Ap will generate a dy- 
namics with a stationary probability distribution, and 
requires 2n + 1 parameters for the symmetric part (the 
elements of Qp, Eq. (HH])), and n{n+l)/2 for the antisym- 
metric part Ap"^^ . If we want the equilibrium distribution 
to be the canonical, one we must enforce the FDT, and 
BpBj is uniquely determined. 



If we aim at a generalized formulation, which allows for 
frequency-dependent thermalization, there are no con- 
straints on the choice of Bp other than the fact that both 
BpBj and the covariance Cp must be positive-definite. 
Clearly, a real, lower-triangular Bp is the most gen- 
eral parameterization of a positive-definite BpB^ , and 

amounts at introducing (n-|-l)(n-|-2)/2 extra parameters. 

(s) 

Together with the assumption that Ap > 0, the con- 
dition BpBp > is sufficient to ensure that the unique 
symmetric Cp which satisfies ([6|) is also positive-definite. 



B. Fitting for canonical sampling 

Armed with such a robust and fairly general parame- 
terization, one only needs to define a merit function to 
be optimized. Again, we first consider the simpler case 
of canonical sampling. Here, we want to obtain a flat re- 
sponse over a wide, physically-relevant frequency range 
i^min , i-^max ) ■ Wc havc choscn the form 



Xi = 



iiogK(wi)r 



(19) 



where w^s are equally spaced on a logarithmic scale over 
the fitted range. If a large value of m is chosen, the 
which yields the lowest efficiency is weighted more, and 
a flat response curve is obtained. We found empirically 
that values of m larger than 10 lead to a proliferation of 
local minima, and hinder efficient optimization. To re- 
solve this, one can use the optimal parameters for m = 2 
as input for further refinement at larger to, until conver- 
gence is achieved. 

This procedure can be modified so as to provide an 
efficient thermostat which can be used in Car-ParrincUo- 
like dynamics. In this case, the GLE has to act as a 
low-pass filter in which only the low ionic frequencies are 
affected, and fast electronic modes are not perturbed. 
To obtain this effect, we compute ([19]) only for the Wi's 
which are smaller than a cutoff frequency ujcp-i and we 
introduce an additional term 



X2 



|max[logK(a;^) - fc(wcp - Wi),0]|^ 



1/m 



(20) 



X2 enforces a steep decrease of k{uj) above (jJcp, with a 
slope A: on a logarithmic scale. Values of k as large as 9 
can be used, which guarantee an abrupt drop in thermal- 
ization efficiency for the fast modes (see Figure [3]). 



C. Non-thermal noise and quantum thermostat 

We now discuss the case in which the thermostat is 
permitted to violate FDT, in order to achieve frequency- 
dependent equilibration. For these applications, one 
must also fit the fluctuations Cpp(a;) and Cgg(w) to some 
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10"' 1 
o) [arb. units] 



FIG. 3: Thermostatting efficiency, as estimated from Eq. ((9|, 
for a colored-noise thermostat optimized for Car-Parrinello 
dynamics. Sampling efficiency is optimized for lj G (10^'^, 1), 
and an abrupt drop in efficiency is enforced for £ (1, 10), 
using the penalty function (|20[l in the fitting. The continuous 
(dark red) curve corresponds to = 9, the dashed (orange) 
curve to = 6 and the dotted (light orange) curve to = 3. 
The K(aj) curve for a white-noise thermostat centered on the 
optimized range is also reported for reference (dotted black 
curve). The three curves correspond to the parameters set 
|cp-9_4-3[ |cp-6_4-3| and ,cp-3_4-3p. 



target function Cpp and Cqq . We shall not treat the general 
case, but rather investigate the example of the quantum 
thermostat (Refi^). The procedure followed provides a 
clear guide for future extensions to different applications. 

In order to reproduce quantum ions effects, one must 
selectively heat high-frequency phonons, for which zero- 
point energy effects are important, without affecting the 
low-frequency modes which behave classically. The re- 
quired frequency dependence of the variance for this case 



is that of a quantum oscillator, i.e. Cppico) 
%^coth 



2k 



Since the low- frequency limit is already enforced by (|2ip , 
we compute (|22p on a set of points equally spaced be- 
tween the maximum frequency LOmax and one half of the 
onset frequency for quantum effects iOq = ksT /h. 



D. Transferability of fitted parameters 

The scheme described in the previous Sections allowed 
us to obtain matrices suitable for all the applications 
discussed in previous works. Furthermore, it provides 
a starting point for obtaining matrices which one might 
deem useful for novel applications. However, the reader 
is advised that the fitting is still far from being a black- 
box procedure. It is thus necessary to experiment with a 
combination of different initial parameters and minimiza- 
tion schemes. We found the downhill simplex methods^^ 
to be particularly effective, but resorted to simulated an- 
nealing when the optimization got stuck in a local mini- 
mum. There is a great deal of arbitrariness in the choice 
of the terms (|19ll22p . and in their weighted combination 
X = ^ WiXi- To make the procedure even more delicate, 
we observe that in high-n cases the parameters tend to 
collapse into "degenerate" minima, where the full dimen- 
sionality of the search space is not exploited. This phe- 
nomenon can be successfully circumvented by enforcing 
an even spacing of the eigenvalues of A over the frequency 
range of interest, and slowly releasing this restraint dur- 
ing the later stages of optimization. 

However, the problems mentioned above have no major 
practical consequences, as the computation of analytical 
estimates is inexpensive and one can afford a great deal of 
trial-and-error during the optimization. Moreover, fitted 
parameters can be reused, since the optimized parame- 
ters can be easily transferred to similar problems because 
of the scaling properties of the dynamics ([5]) . 

In fact, one can see that if the drift and covariance 
matrices (Ap,Cp) lead to the efficiency curves k(uj) 
and fluctuations Cpp{uj), the scaled matrices (a Ap, /? Cp) 



^ The w — > 0, classical limit can be proved to -i, • , i / _i f ' j ii, a j. i- a t -i \ 
T '. fir yield K.{a uj), and the fluctuations pcpp(a ^w). 



correspond to two conditions on the elements of the free- 
particle covariance matrix Cp; namely, Cpp = fesT and 



apA- 



Cp = 0. One could enforce such constraints ex- 



actly, by considering the entries of Cp as independent fit- 
ting parameters, and obtaining the diffusion matrix from 
Eq. We found however that this choice makes it dif- 
ficult to obtain a positive-definite BpBj, and that the 
fitting becomes more complex and inefficient. 

As an alternative, we decided to enforce the low- 
frequency limit with an appropriate penalty function. 



X3 = (cpp/kBT - 1)2 + (ajA-^Cp/fcBT) 



(21) 



to be optimized together with the sampling efficiency 
and a term which measures how well the finite-frequency 
fiuctuations were fitted: 



X4 = 



E 



log- 



qq 



+ 



log Cpp (w,;) 



1/m 



(22) 



This means that if Ap is optimized for sampling over 
the range {ujrninT^max), aAp will be optimal over 
(a ujmim ct ujmax) ■ Wc also remark that if (Ap,Cp) are 
fitted to the quantum harmonic oscillator fluctuations at 
temperature T, (a Ap, a Cp) will be suitable for tempera- 
ture a T. Care must be taken in this case to ensure that 
the scaled frequency range still encompasses the whole 
vibrational spectrum of the system being studied. 



III. UNDERSTANDING THE QUANTUM 
THERMOSTAT 

As discussed in Ref<^, one must pay a great deal of at- 
tention when using a "quantum thermostat" , because en- 
ergy is transferred between modes of different frequency, 
as a consequence of the anharmonic coupling. This is 
reminiscent of zero-point energy (ZPE) leakage which 
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plagues semiclassical approaches to the computation of 
nuclear quantum effects^i^. In the cases we explored 
so far, empirical evidence suggests that quasi-harmonic 
solids, can be treated with good accuracy down to tem- 
peratures as low as 10% of the Debye temperature Qd- 
Clearly, the ultimate test to assess of the accuracy of 
the method is a comparison with path-integral calcula- 
tions, to be performed on a similar but computationally 
cheaper model, such as a smaller-size box or a simpler 
force field. 

One would like however to obtain some qualitative 
measure of the quality of the fit, and to gauge the trans- 
ferability of a given set of parameters. To this end, we 
first state a couple of empirical rules, and then validate 
them on two fairly different real systems. A first obser- 
vation is that it is useless to push the fitting of the fluc- 
tuations Cpp{uj) and Cqq{io) to very high accuracy, if this 
comes at the expense of the coupling efficiency. In fact, 
we would be trading a small, controlled fitting error with 
a possibly larger, uncontrollable and system-dependent 
error stemming from anharmonicity. Secondly, we ob- 
served that in order to contrast more effectively the flow 
of energy between different phonons, one should try to re- 
duce the correlation time of the kinetic energy tk, rather 
than focus solely on the terms which are better suited 
to measure sampling efficiency. In fact, a low tx(w) cor- 
responds to a slightly overdamped regime, where sam- 
pling efficiency is sub-optimal, but ZPE is enforced more 
tightly 

To demonstrate these concepts in a real system, we 
performed some calculations with a Tersoff model of di- 
amond at a temperature T = 200 K. At this low tem- 
perature, slightly below O.IO^), quantum effects are very 
strong, and we therefore expect to have problems main- 
taining the large difference in temperature between the 
stiff and soft phonons. Using a very harmonic system 
such as diamond is particularly useful, since one can mon- 
itor directly the efficiency of the thermostat by project- 
ing the atomic velocities on a selection of normal modes. 
Hence, a projected kinetic temperature T'{uj) can be 
computed, and its value checked against the predictions 
in the harmonic limit, in the same spirit as in Rcfj^. In 
Figure [4] we report the results with a matrix fitted tak- 
ing into account only the terms (pij) and (P^. Even in 
a harmonic system such as diamond there are major er- 
rors due to ZPE leakage from the high-frequency to the 
low- frequency modes, which the thermostat compensates 
only partially. These poor results should be compared 
with those of Figure [5] Here, we have also introduced 
in the fit a term analogous to to reduce the value 
of r/f (w). The projected kinetic temperature now agrees 
almost perfectly with the analytical predictions Cpp(w) 
for most of the modes. The only ones displaying signif- 
icant deviations are the faster ones, for whom the value 
of Tif(a;) is slightly larger. The Cpp(w) curve deviates 
by nearly 10% from the exact, quantum-mechanical ex- 
pectation value. However, thanks to the more efficient 
coupling, the errors due to anharmonicities are better 
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FIG. 4: (a): oj-dependence of the kinetic energy correlation 
time Tk{uj) (light, dotted line) and of the ratio of the fit- 
ted fluctuations Cpp{uj) (dashed line) and of Lo'^Cqq{uj) (full 
line) with the exact, quantum-mechanical target function, 
(b): normal-mode-projected kinetic temperature for a few, 
selected phonons. The dashed line is the value expected 
from the fitting Cpp(aj), while the full line is the exact, 
quantum-mechanical expectation value for a harmonic oscil- 
lator. Calculations have been performed with the parameters 
,qt-20_6_BADf^. 
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FIG. 5: (a): cj-dependence of the kinetic energy correla- 
tion time T_f(-(tj) (light, dotted line) and of the ratio of 
the fitted fluctuations Cpp{Lj) (dashed line) and of ui'^Cqq{u)) 
(full line) with the exact, quantum-mechanical target func- 
tion, (b): normal- mode-projected kinetic temperature for 
a few, selected phonons. The dashed line is the value ex- 
pected from the fitting Cpp{ijj), while the full line is the exact, 
quantum-mechanical expectation value for a harmonic oscil- 
lator. Calculations have been performed with the parameters 
qt-20_6J^. 
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FIG. 6: Radial distribution function as computed from 
fully-converged path-integral calculations^^ (black, dotted 
line), and from a quantum-thermostat MD trajectory for a 
Lennard- Jones model of solid neon at T = 20 K. Distances 
are in reduced units. Full line corresponds to the parameters 
set |qt-20_6 | (cfr. Figure [Sjl , and lighter, dashed line to the 
set |qt-20_6_BAD| (cfr. Figure |4j. 

compensated, and in actuality, the overall error is much 
smaller than for the parameters presented in Figure ID 

To test whether these prescriptions work for less har- 
monic problems, we now turn to a completely different 
system; namely, the structural properties of solid neon 
at 20 K. At variance with diamond, quantum-ions effects 
are less pronounced, but the system is close to its melt- 
ing temperature, and it is significantly anharmonic. As 
shown in Figure |6l the agreement between our results 
and those of accurate path-integral calculations^ is al- 
most perfect if the parameters of Figure [5] are used. As 
expected, large errors are present if|q t-20_6_BAD| is used. 
Further improvements on the fitting strategy, and the ap- 
plication to strongly anharmonic systems is currently be- 
ing investigated, and will be the subject of further work. 



IV. CONCLUSIONS 



from one system to another. With this in mind we have 
provided an extensive library of optimized parametersii, 
which makes fitting unnecessary for most applications. 

We also comment on practical issues concerning the 
implementation of the generalized-Langevin thermostat 
in a molecular-dynamics program and its use in appli- 
cations. In particular, we discuss in detail how one can 
use colored noise to model nuclear quantum effects^. We 
provide some empirical rules to guide the fitting in this 
difficult case, and we demonstrate that a normal-mode 
analysis in a quasi-harmonic system is a valuable tool 
for assessing the quality of a set of parameters. We be- 
lieve that further investigation will find many other appli- 
cations for colored-noise in molecular-dynamics, and in 
computer simulations of molecular systems in general. As 
an example, we are currently investigating using a zero- 
temperature, optimal-sampling GLE thermostat in order 
to perform structural optimization. On similar lines, and 
taking inspiration from "quantum annealing" '^^i^^ , one 
can envisage using frequency-dependent thermalization 
to improve the performance of simulated annealing. 
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In this paper we have discussed in detail the use of 
colored-noise dynamics, based on Ornstein-Uhlenbeck 
processes, as a tool for performing molecular dynam- 
ics. Applications range from enhanced sampling, which 
we demonstrate in the harmonic limit and which will 
be applied to real systems in forthcoming publications, 
to thermostats for adiabatically separated problems and 
frequency-dependent thermalization. 

Our idea exploits the linear nature of the OU stochas- 
tic differential equations, which allows one to use the 
one-dimensional harmonic oscillator as a simple but 
physically- motivated reference model. On the basis of the 
analytical prediction obtained in that case, we describe a 
recipe for fitting the thermostat parameters so as to ob- 
tain the desired response properties in real systems. The 
procedure is not simple, and we are considering different 
approaches to make it more robust and effective. Fortu- 
nately however, fitted matrices can be easily transferred 



Appendix A: Memory kernels for the 
non-Markovian formulation 



The connection between the Markovian ([2]) and non- 
Markovian (O formulations of the colored- noise Langevin 
equation can be understood using techniques similar to 
those adopted in Mori-Zwanzig theorjj^iii. Let us first 
consider a very general, multidimensional OU process, 
where we single out some degrees of freedom (y) that we 
wish to integrate out, leaving only the variables marked 
as x. 

x_\ _ / A^JA^\ / x\ f_^xi_ 

Assuming that the dynamics has finite memory, one 
can safely take y(— oo) = 0, and the ansatz 



(Al) 



y{t) = f e-(*-*')^- [-AyMt')+'Byim]dt'. 

J — OO 



(A2) 
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Substituting into (jAl[) . one sees that y can be eliminated 
from the dynamics of x, and arrives at 



±{t) = - I K(t - t')x(t')dt' + C{t) 

J —OO 

K{t) ^2A,,S{t) ~ A^yc-^'^yy Ay, (t > 0) (A3) 
C(t) =B,5^(t) - [ A,,e-(*-*')^-B,5|(<'). 

^ — OO 



One can see that (|A3P are invariant under any orthogo- 
nal transformation of the y dynamical variables, meaning 
that such a transformation leaves the dynamics of the x's 
unchanged. 

The colored noise is better described in terms of its 
time-correlation function, H(t) = (C(i)C(O)"^)- Let 
us first introduce the symmetric matrix D = BB"^, 
whose parts we shall label using the same scheme 
used for A in Eq. (jAip . We shall also need Zyy = 

Io° e~'^'"'*Dyj,e~^!'!'*dt. With these definitions in mind, 
one finds 

H(i) = J(t)D,, + A^^e-*^- [Zyy^ly - ^yx] {t > 0) . 

(A4) 

Note that the value of H(t) for i < is determined by 
the constraint H(— i) = H(t)-^; the value of K(<) in- 
stead, is irrelevant for negative times: we will assume 
K(— i) = K(t)-'" to hold, since this will simplify some 
algebra below. 

Let's now switch to the case of the free-particle coun- 
terpart of Eqs. ([2]), which is relevant to the memory func- 
tions entering Eqs. ([7]). Here, we want to integrate away 
all the s degrees of freedom, retaining only the momen- 
tum p. Hence, we can transform Eqs. (|A3[) and (|A4[) to 
the less cumbersome form 



-|*|A, 



K{t) =2appS{t) 
Hit) ^dpp6{t) - aje-l*l^ [Za^ - dp] 



(A5) 



This compact notation hides certain relevant property 
of the memory kernels, which are more apparent when 
the kernels are written in their Fourier representation. If 
Dp = BpBj is transformed according to Eq. K{uj) 
and H{uj) read 



K{uj) ~2app — 2a, 



H{uj) =K{lj) ( Cpp 
1 



f A2+C^2-P 



(A6) 



20.2 



> A 2 I ,,,2^P 



l + a. 



V A2 



,2"P 



It is seen that the memory functions (hence the dynam- 
ical trajectory) are independent of the value of C, the 
covariance of the fictitious degrees of freedom. More- 
over, a sufficient condition for the FDT to hold is readily 
found. By setting Cpp = kgT and Cp = 0. one obtains 
H{ijj) ~ kgTKiuj), which is precisely the FDT for a non- 
Markovian Langevin equation. Since the value of C is 



irrelevant we can take Cp = ksT, which simplifies the 
algebra and leads to numerically-stable trajectories. 



Appendix B: Covariance matrix and correlation 
times for the harmonic oscillator 



Given A and C matrices (the drift term and the static 
covariance for a generic OU process), one can find the dif- 
fusion matrix B by an expression analogous to Eq. ©. 
The same relation can be used to obtain the elements of 
C given the drift and diffusion matrices, by solving the 
linear system. However, the covariance matrix can be 
computed more efficiently by finding the eigendecompo- 
sition of A = O diag(Q!i) 0~^, and computing 



O-iBB^O 



O 



kl 



o-k + ai 



(Bl) 



Now, let X be the vector describing the trajectory 
of the OU process. In order to compute th or ry 
(Eq. ([9|)) one needs time correlation functions of the 
form {xi{t)xj{t)xk{0)xi(0)) . The corresponding, non- 
normalized integrals 



Tijkl 



{xi{t)xj{t)xk{0)xi{0)) - {x,Xj) {xkXi)]dt 

(B2) 



can be computed in terms of the tensorial quantity 

0„n [O-lC" 



ijkl 



■il J"- ink 



(B3) 



as Tijkl = \ {Xijki + Xijik + Xkiij + Xikij). For example 
~ if we consider the full OU process in the harmonic case 
- one computes 



TH = 



iqqqq 



qqpp 



' pppp 



TV 



Iqqqq 



(B4) 



PP n 
where we use an obvious notation for the indices in t,, 



ijkl • 



Appendix C: A comparison with Nose-Hoover 
Chains 



The most widespread techniques for canonical sam- 
pling in MD are probably white-noise Langevin and Nose- 
Hoover chains (NHC). White- noise Langevin can be con- 
sidered as a limit case of the thcrmostatting method we 
describe in this work, but NHC is based on a redically 
different philosophy. It is therefore worth performing a 
brief comparison between the latter and the GLE ther- 
mostat. 

In the "massive" version of the NH thermostaij^^— , 
each component of the physical momentum is coupled 
to an additional degree of freedom with a fictitious mass 
Q, by means of a second-order equation of motion. The 
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resulting dynamics ensures that the physically-relevant 
degrees of freedom will sample the correct, constant- 
temperature ensemble, with the advantage of having de- 
terministic equations of motion, and a well-defined con- 
served quantity. However, in the harmonic case, trajecto- 
ries are poorly ergodic. This problem can be addressed by 
coupling the fictitious momentum to a second bath vari- 
able with a similar equation of motion. By repeating this 
process further a "Nose-Hoover chain" can be formed, 
which ensures that the dynamics is sufficiently chaotic to 
achieve efficient samplingi^i^. The drawback of this ap- 
proach is that the thermostat equations are second-order 
in momenta. It is therefore difficult to obtain analyti- 
cal predictions for the properties of the dynamics, and 
the integration of the additional degrees of freedom must 
be performed with a multiple time-step approach, which 
makes the thermostat more expensive. 

To examine the performances of NHC and GLE, one 
could envisage comparing the sampling efficiency as de- 
fined by the correlation times (0). Obtaining such esti- 
mates is not straightforward, not only because the the 
harmonic case cannot be treated analytically, but also 
because in the multidimensional case the properties of 
the trajectory will not be invariant under an orthogonal 
transformation of coordinates, as discussed in Section HI 
The simplest model we can conceive for comparing NHC 
and GLE is therefore a two-dimensional harmonic oscil- 



lator, with different vibrational frequencies on the two 
normal modes and adjustable relative orientations of the 
eigenvectors with respect to the thermostatted coordi- 
nates. 

The resulting Ty is reported in Figure [T] in the highly 
anisotropic cases, the efficiency of the NH chains depends 
dramatically on the orientation of the axes, while for 
well-conditioned problems is almost constant. The lin- 
ear stochastic thermostat, on the other hand, has a pre- 
dictable response, which is completely independent on 
orthogonal transforms of the coordinates. In the one- 
dimensional case - or when eigenvectors are perfectly 
aligned with axes - NH chains are very efficient for all 

modes with frequency w < \J^^-- One should how- 
ever consider that, in the absence of an exact propaga- 
tor, choosing a small Q implies that integration of the 
trajectory for the chains will become more expensive. 

Obviously, such a simple toy model does not give quan- 
titative information on the behavior in real-life cases, 
where modes of different frequencies coexist with an- 
harmonicity and diffusive behavior. However, it demon- 
strates that the colored-noise Langevin thermostat per- 
forms almost as well as the axis-aligned NH chains. Fur- 
thermore, unlike the NHC, there are no unpredictable 
failures for anisotropic potentials. 
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FIG. 7: Correlation time for the potential energy of a 2-D 
harmonic oscillator, as a function of the angle between the 
eigenmodes and the cartesian axes, tv is computed for dif- 
ferent values of the condition number uimax/ij^min, from bot- 
tom to top 10, 31.6 and 100. Thin lines serve as an aid for 
the eye, connecting the results obtained in the three cases 
using a massive NH chains thermostat, with four additional 
degrees of freedom and Q — ksT /lu^^^. Error bars are also 
shown for individual data points. Thick lines correspond to 
the (constant) result predicted for a CLE thermostat, using 
respectively the thermostat parameters kv_2- 1 , centered on 
Q.i2ujmax, kv_4-2, centered on Q.l^ujrnax, and kv_4-2 cen- 
tered on Q.lu)rnax- The values obtained in actual GLE simu- 
lations agree with the predictions within the statistical error- 
bar, and are not reported. 



